﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
//using System.Drawing;

using Qiang.Topology.Geometries;


namespace Qiang.QMath
{
    public class QMath
    {

        //地理信息系统算法基础。
        // 我还没看懂的有：第三章（空间数据变换算法）、第四章的栅格转矢量、



        //设矢量P=（X1,Y1），Q=（X2,Y2），则矢量叉积定义为由（0，0）(原点其实可自定义)、P、Q和PQ所组成的平行四边形的带符号的面积。
        //即：P*Q=X1*Y2-X2*Y1,其结果是一个标量（只有大小没有方向）。
        //注意：是有先后顺序的
        //若P*Q>0，则P在Q的顺时针方向，=0则共线，<0则在逆时针方向


        //判断点Q是否在直线P1P2上，则直接考虑矢量P1Q与P1P2是否共线即可。在这里可认为P1点为原点。矢量P1Q为（Q.X-P1.X，Q.Y-P1.Y）,同理。
        //再判断点Q是不否在线段P1P2上，直接考虑点是Q是否在P1P2为对角线的矩形区域内


        /// <summary>
        /// 矢量法判断点是否在线上。也可以此方法判断点是在线（矢量线）的哪个方向（顺时针or逆时针）
        /// </summary>
        /// <param name="pt"></param>
        /// <param name="line"></param>
        public static bool PtInLine(Coordinate pt, Coordinate[] line)
        {
            //矢量1
            Coordinate v1 = new Coordinate(pt.X - line[0].X, pt.Y - line[0].Y);
            //矢量2
            Coordinate v2 = new Coordinate(line[1].X - line[0].X, line[1].Y - line[0].Y);
            //矢量相乘判断是否共线
            double b = 0;
            b = v1.X * v2.Y - v2.X * v1.Y; //矢量叉积判断是否共线
            //if (b == 0)  //要考虑误差
            if (Math.Abs(b) < 0.000001) //等于0（在误差范围内）共线
            {
                //判断是否可能在延长线上
                if (pt.X <= Math.Max(line[0].X, line[1].X) && pt.X >= Math.Min(line[0].X, line[1].X) &&
                    pt.Y <= Math.Max(line[0].Y, line[1].Y) && pt.Y >= Math.Min(line[0].Y, line[1].Y))
                {
                    return true;
                }
            }
            return false;

        }

        /// <summary>
        /// 判断点在线的哪个方向，等于0说明在线上或线的延长线上，大于0在线的右手边（顺时针方向），小于0在左手边
        /// </summary>
        /// <param name="pt"></param>
        /// <param name="line"></param>
        /// <returns></returns>
        public static double PtAtLine(Coordinate pt, Coordinate[] line)
        {
            //矢量1
            Coordinate v1 = new Coordinate(pt.X - line[0].X, pt.Y - line[0].Y);
            //矢量2
            Coordinate v2 = new Coordinate(line[1].X - line[0].X, line[1].Y - line[0].Y);
            //矢量相乘判断是否共线
            double b = 0;
            b = v1.X * v2.Y - v2.X * v1.Y; //矢量叉积判断是否共线

            if (Math.Abs(b) < 0.000001) b = 0;  //误差判断

            return b;
        }


        /// <summary>
        /// 判断两条线段AB与CD是否相交。
        /// 矢量叉积法，可判断是相交、相切、相离。但是不知道怎么求交点坐标。
        /// 此方法精确，而且没有除法算法。
        /// </summary>
        /// <param name="AB"></param>
        /// <param name="CD"></param>
        /// <returns></returns>
        public static bool LineIntersect(Coordinate[] AB, Coordinate[] CD)
        {
            //先判断两条线做为对角线所在的矩形区域是否相交
            Coordinate[] extent1 = new Coordinate[2];
            extent1[0] = new Coordinate(Math.Min(AB[0].X, AB[1].X), Math.Min(AB[0].Y, AB[1].Y));
            extent1[1] = new Coordinate(Math.Max(AB[0].X, AB[1].X), Math.Max(AB[0].Y, AB[1].Y));

            Coordinate[] extent2 = new Coordinate[2];
            extent2[0] = new Coordinate(Math.Min(CD[0].X, CD[1].X), Math.Min(CD[0].Y, CD[1].Y));
            extent2[1] = new Coordinate(Math.Max(CD[0].X, CD[1].X), Math.Max(CD[0].Y, CD[1].Y));

            //先快速判断
            if (extent1[0].X > extent2[1].X || extent1[0].Y > extent2[1].Y || extent1[1].X < extent2[0].X || extent1[1].Y < extent2[0].Y)
            {
                return false;
            }
            //通过了快速判断再进行细致判断
            //如果线AB与CD相交，则C、D两点必定一个在A的顺时针方向，一个在A的逆时针方向，
            //根据叉积的性质：矢量AC和AB的叉积 与 AD与AB的叉积必定一个是正数一个是负数。两个叉积相乘（叉积是标量）必定为负数。
            //此处设A为原点来构造矢量……

            //矢量AC
            Coordinate vectorAC = new Coordinate(CD[0].X - AB[0].X, CD[0].Y - AB[0].Y);
           //矢量AD
            Coordinate vectorAD = new Coordinate(CD[1].X - AB[0].X, CD[1].Y - AB[0].Y);
            //矢量AB
            Coordinate vectorAB = new Coordinate(AB[1].X - AB[0].X, AB[1].Y - AB[0].Y);
            //AC与AB的叉积
            double CAB = vectorAC.X * vectorAB.Y - vectorAC.Y * vectorAB.X;
            //AB与AD的叉积  //第二个叉积反着来，矢量AB在前，CD在后，这样两个叉积相乘为正才是相交
            double BAD = vectorAB.X * vectorAD.Y - vectorAB.Y * vectorAD.X;

            //if (CAB == 0 && BAD == 0) return 重合;
            //if (CAB ==0 && BAD!=0) return C在线AB上，有可能各A点或B点重合。
            //if(BAD ==0&& CAB !=0) return D在线AB上，有可能各A点或B点重合。

            double resultCD = CAB * BAD;  //result大于0只能证明CD跨立AB（有可能A在线CD上，或B在CD上），还得证明AB跨立CD

            if (resultCD > 0)
            {
                //再以C点为原点证明AB跨立CD
                Coordinate vectorCA = new Coordinate(-vectorAC.X, -vectorAC.Y);
                Coordinate vectorCB = new Coordinate(AB[1].X - CD[0].X, AB[1].Y - CD[0].Y);
                Coordinate vectorCD = new Coordinate(CD[1].X - CD[0].X, CD[1].Y - CD[0].Y);

                double ACD = vectorCA.X * vectorCD.Y - vectorCA.Y * vectorCD.X;
                double DCB = vectorCD.X * vectorCB.Y - vectorCD.Y * vectorCB.X;

                //if (ACD == 0 && DCB == 0) return 重合;
                //if (ACD == 0 && DCB != 0) return A点在线CD上，有可能和A、B重合;
                //if (DCB == 0 && ACD != 0) return B点在线CD上，有可能和A、B重合;

                double resultAB = ACD * DCB;
                if (resultAB > 0) return true; //AB也跨立CD

            }
            

            return false;
        }


        /// <summary>
        /// 射线法判断点是否在多边形内。用到了两条线段是否相交的函数方法。
        /// 这个方法比以前写的那个更精确更稳定更快捷，只需做乘法和减法。
        /// 还有一个转角法，暂时不会写……
        /// </summary>
        /// <param name="pt"></param>
        /// <param name="polygon"></param>
        /// <returns></returns>
        public static bool PtInPolygon(Coordinate pt, Coordinate[] polygon)
        {
            Coordinate minPt = new Coordinate(double.MaxValue, double.MaxValue);
            Coordinate maxPt = new Coordinate(double.MinValue, double.MinValue);

            foreach (Coordinate vertex in polygon)
            {
                //快速判断点是否与多边形的顶点重合
                if (pt.X == vertex.X && pt.Y == vertex.Y)
                {
                    return false;
                }
                //计算多边形的范围（boundingBox）
                if (minPt.X > vertex.X)
                {
                    minPt.X = vertex.X;
                }
                if (minPt.Y > vertex.Y)
                {
                    minPt.Y = vertex.Y;
                }
                if (maxPt.X < vertex.X)
                {
                    maxPt.X = vertex.X;
                }
                if (maxPt.Y < vertex.Y)
                {
                    maxPt.Y = vertex.Y;
                }
            }
            //快速判断点是否在多边形的范围内
            if (pt.X < minPt.X || pt.Y < minPt.Y || pt.X > maxPt.X || pt.Y > maxPt.Y)
            {
                return false;
            }

            //===========================以上是快速判断===================================//
            //==========================通过了快速判断，射线法精确判断====================//
            //构造向右的水平射线线段
            Coordinate[] rayLine = new Coordinate[2];
            rayLine[0] = new Coordinate(pt.X, pt.Y);
            rayLine[1] = new Coordinate(maxPt.X+10, pt.Y);

            int num = 0;

            for (int i = 0; i < polygon.Length - 1; i++)
            {
                Coordinate[] line = new Coordinate[2];
                line[0] = new Coordinate(polygon[i].X, polygon[i].Y);
                line[1] = new Coordinate(polygon[i + 1].X, polygon[i + 1].Y);

                if (line[0].Y == line[1].Y) continue; //水平的线不参与判断
                //如果射线穿过多边形某条线的端点，按以下规则进行处理：
                //方向向上的边包括其开始点，不考虑终止点
                //方向向下的边包括其终止点，不考虑开始点
                if (line[1].Y > line[0].Y)
                {
                    if (line[0].X > pt.X)
                    {
                        if (line[0].Y == pt.Y)
                        {
                            num++;
                            continue;
                        }
                    }
                    if (line[1].Y == pt.Y) continue;
                }
                else
                {
                    if (line[0].Y == pt.Y) continue;
                    if (line[1].X > pt.X)
                    {
                        if (line[1].Y == pt.Y)
                        {
                            num++;
                            continue;
                        }
                    }
                }
                //特殊情况判断完后，判断普通情况两条线是否相交
                bool b = LineIntersect(rayLine, line);
                if (b == true) num++;
            }
            int n = 0;
            Math.DivRem(num, 2, out n);
            //交点个数是奇数则在多边形内
            if (n == 1) return true;


            return false;
        }

        

        // 矢量线的栅格化有3种方法
        // 1 八方向栅格化。
        //   根据矢量的倾角情况，在每行每列上，只有一个像元被涂黑（栅格化）。
        //   假如行数大于列数（即矢量线两点Y的差值大于两点X的差值），则逐行求出本行中心线与线的交点，交点所在像元涂黑。反之列大于行。
        // 2 全路径栅格化。
        //   即线穿过哪个网格，哪个网格就被涂黑。
        //   书上的方法没看懂。不过可以用线与网格相交的方法来判断啊……
        // 3 恒密度栅格化。（书上没讲。。。）

        /// <summary>
        /// 矢量到栅格（八方向栅格化）。伪代码。。。
        /// </summary>
        /// <param name="line"></param>
        public static void VectorToRaster(Coordinate[] line)
        {
            //栅格图像的大小、行列数。行列数大者默认为255
            int row = 0;
            int col = 0;
            double stepX = 1;
            double stepY = 1;
            int[,] data;  //二维栅格数据，自己虚构（有数据存储1，无数据存储无效值（默认为0）。。。

            //先判断行列数
            if (line[1].X - line[0].X < line[1].Y - line[0].Y)
            {
                row = 255;                

                if (Math.Abs(line[1].Y - line[0].Y) < 0.1)  //得用误差。线是水平的情况
                {
                    row = 1;
                    col = 1;
                }
                else
                {
                    col = (int)Math.Floor((line[1].X - line[0].X) * 255 / (line[1].Y - line[0].Y)) + 1;  
                }
                //具有方向的步长。(有可能为0，不加判断了，按一般情况处理吧
                stepX = (line[1].X - line[0].X) / col;
                stepY = (line[1].Y - line[0].Y) / row;

                //初始化数据
                data = new int[row, col];                
                for (int i = 0; i < row; i++)
                {
                    for (int j = 0; j < col; j++)
                    {
                        data[i, j] = 0;
                    }
                }
                //计算有效数据 
                double a = (line[1].Y - line[0].Y) / (line[1].X - line[0].X);
                double b = line[0].Y - a * line[0].X;

                for (int i = 0; i < row; i++)
                {
                    double Y = stepY / 2 + i * stepY;

                    if (line[0].X == line[1].X)
                    {

                    }
                    else
                    {
                        double X = (Y-b)/a;     //a是否为0.都要加判断的……
                        int j = (int)Math.Floor(X / stepX);
                        data[i, j] = 1;
                    }

                }

            }
            else
            { 
                //道理同上。。。。
            }
        }


        //矢量面的栅格化（即矢量多边形添充。）
        //栅格转矢量还没看懂。。。。。。。

        //===========================================================================================//

        //矢量线的压缩算法
        // 1 道格拉斯-普克法，类似垂距法，下面有代码
        // 2 光栏法。书上说这个算法运算量少，算法严密，但我感觉这个方法太繁锁，代码量大，还容易出错



        //点到直线的距离公式： 其实有两种方法
        // 1 点（xi,yi）到直线Ax+By+C=0的距离为： d=|Axi+Byi+C|. 这个方法还要求直线的方程系数，A B C。这种方法就是需要判断直线的情况、水平、垂直
        // 2 先求出点P到线两端点L1、L2的距离 d1 d2,和线本身的长度L（直接求两点距离）
        //   再求出两端点L1、L2分别到垂点R的距离d3 d4。这个公式为d3=（d1^2+L^2-d2^2）/(2*L)  d4=(d2^2+L^2-d1^2)/(2*L).这个公式可用三角函数证明.(这两个值是有正负的，关系到夹角的关系，锐角为正，钝角为负)
        //   最后用上面的结果求点到直线的距离d=math.sqrt(L1^2-d3^2),或d=math.sqrt(L2^2-d4^2)

        

        /// <summary>
        /// 矢量线的压缩算法.道格拉斯-普克法。这个结果是反序的
        /// </summary>        
        public static void YaSuoLine(Coordinate[] line,double limitDistance,ref List <Coordinate> result)
        {
            if (line.Length < 2) return;

            //必要的条件判断
            if (limitDistance == 0)
            {
                result.Insert(0, line[0]);
                return;
            }

            result.Insert(0, line[line.Length - 1]);

            //递归结束条件
            if (line.Length == 2)
            {
                result.Insert(0, line[0]);
                return;
            }
                       

            //存储有效的点，再对有效的点用同样的方法进行筛选
            List<Coordinate> temLine = new List<Coordinate>();
            temLine.Add(line[0]);

            double A, B, C;
            if (line[0].X != line[line.Length - 1].X)
            {
                //第一个点到最后一个点的直线方程 Ax+By+C=0
                A = (line[line.Length - 1].Y - line[0].Y) / (line[line.Length - 1].X - line[0].X);
                B = -1;
                C = line[0].Y - A * line[0].X;

                //遍历中间的点，判断中间点到首尾连线的距离进行筛选
                for (int i = 1; i < line.Length - 1; i++)
                {
                    double distance = Math.Abs(A * line[i].X + B * line[i].Y + C);  //点到直线的距离公式
                    if (distance >= limitDistance) //如果大于有效距离，则此点是有效的
                    {
                        temLine.Add(line[i]);
                    }
                }

                //递归调用
                YaSuoLine(temLine.ToArray(), limitDistance, ref result);

            }
            else
            { 
                //原理同上
            }
            
        }

        /// <summary>
        /// 光栏法,矢量线的压缩。暴露给外部
        /// </summary>
        /// <param name="line"></param>
        /// <param name="distance"></param>
        /// <param name="result"></param>
        public static void CompressLine(Coordinate[] line, double distance, ref List<Coordinate> result)
        {
            Coordinate[] lastPt = new Coordinate[2];
            lastPt[0].X = float.NaN;
            lastPt[0].Y = float.NaN;
            lastPt[1].X = float.NaN;
            lastPt[1].Y = float.NaN;
            //lastPt先自己定义
            CompressLine(line, distance, lastPt, ref result);
        }



        /// <summary>
        /// 光栏法 ，矢量线的压缩。内部递归算法，供上面那个方法调用 
        /// </summary>
        /// <param name="line"></param>
        /// <param name="distance">有效距离</param>
        /// <param name="result"></param>
        private static void CompressLine(Coordinate[] line, double distance, Coordinate[] lastPt, ref List<Coordinate> result)
        {
            if (line.Length < 2) return;

            //如果第三个点在前两个点构成的扇形范围内，则会删除第二个点递归调用，以防止重复添加第一个点，这里加判断          
            
            if (result.Contains(line[0]) == false)  
            {
                result.Add(line[0]);
            }
            //如果上面这个判断不行就用另一个判断 
            //if (result.Count == 0 || (result[result.Count - 1].X != line[0].X || result[result.Count - 1].Y != line[0].Y))
            //{
            //    result.Add(line[0]);
            //}
            

            //必要的条件判断
            if (distance == 0)
            {
                result.Add(line[line.Length - 1]);
                return;
            }

            //递归结束条件
            if (line.Length == 2)
            {
                result.Add(line[line.Length - 1]);
                return;
            }

           //下面这个求pt，还有改进的空间，可以不用传pt，直接传参 矢量叉积就行


            double A, B, C;

            Coordinate[] pt = new Coordinate[2];
            //求有效距离与线组成的扇形的两点
            if (line[0].X != line[1].X)
            {
                A = (line[1].Y - line[0].Y) / (line[1].X - line[0].X);
                B = -1;
                C = line[0].Y - line[0].X * A;

               
                if (line[0].Y != line[1].Y)
                {
                    pt[0].X = (float)((distance / 2 - B * line[1].X / A - B * line[1].Y - C) / (A - B / A));
                    pt[0].Y = (float)(line[1].Y - (pt[0].X - line[1].X) / A);
                    pt[1].X = (float)((-distance / 2 - B * line[1].X / A - B * line[1].Y - C) / (A - B / A));
                    pt[1].Y = (float)(line[1].Y - (pt[1].X - line[1].X) / A);
                }
                else
                {
                    pt[0].X = line[1].X;
                    pt[0].Y = (float)(line[1].Y - distance / 2);
                    pt[1].X = line[1].X;
                    pt[1].Y = (float)(line[1].Y + distance / 2);
                }

            }
            else
            {
                pt[0].X =(float)( line[1].X - distance / 2);
                pt[0].Y = line[1].Y;
                pt[1].X = (float)(line[1].X + distance / 2);
                pt[1].Y = line[1].Y;

            }


            //与上一次的范围进行比较
            if (lastPt[0].X != float.NaN && lastPt[0].Y != float.NaN && lastPt[1].X != float.NaN && lastPt[1].Y != float.NaN)
            {
                //矢量叉积判断这四个点在线的哪一侧，哪两个点与线构成的角度更小
                double a = (lastPt[0].X - line[0].X) * (line[1].Y - line[0].Y) - (lastPt[0].Y - line[0].Y) * (line[1].X - line[0].X);
                double b = (lastPt[1].X - line[0].X) * (line[1].Y - line[0].Y) - (lastPt[1].Y - line[0].Y) * (line[1].X - line[0].X);
                double c = (pt[0].X - line[0].X) * (line[1].Y - line[0].Y) - (pt[0].Y - line[0].Y) * (line[1].X - line[0].X);
                double d = (pt[1].X - line[0].X) * (line[1].Y - line[0].Y) - (pt[1].Y - line[0].Y) * (line[1].X - line[0].X);

                if (a < b)
                {
                    if (c < d)
                    {
                        if (a > c)
                        {
                            pt[0] = lastPt[0];
                        }
                        if (b < d)
                        {
                            pt[1] = lastPt[1];
                        }
                    }
                    else
                    {
                        if (a > d)
                        {
                            pt[1] = lastPt[0];
                        }
                        if (b < c)
                        {
                            pt[0] = lastPt[1];
                        }
                    }
                }
                else   //a>b
                {
                    if (c > d)
                    {
                        if (a < c)
                        {
                            pt[0] = lastPt[0];
                        }
                        if (b > d)
                        {
                            pt[1] = lastPt[1];
                        }
                    }
                    else
                    {
                        if (a < d)
                        {
                            pt[1] = lastPt[0];
                        }
                        if (b > c)
                        {
                            pt[0] = lastPt[1];
                        }
                    }
                }

                //矢量叉积判断第三个点在线的哪一侧，并用值来判断是否在扇形范围内
                double isIn = (line[2].X - line[0].X) * (line[1].Y - line[0].Y) - (line[2].Y - line[0].Y) * (line[1].X - line[0].X);

                List<Coordinate> ptlst = new List<Coordinate>();

                if (isIn > Math.Min(c, d) && isIn < Math.Max(c, d))  //第三个点前两个点与给定有效距离所形成的扇形之间，则删除第二个点，以第一个点为起点重新递归调用
                {
                    ptlst.AddRange(line);
                    ptlst.RemoveAt(1);  //移除第二个点
                    CompressLine(ptlst.ToArray(), distance, pt, ref result);  //pt做为上一次留下的

                }
                else  //第三个点不在范围内，删除第一个点，从第二个点开始递归进行下一次判断
                {
                    ptlst.AddRange(line);
                    ptlst.RemoveAt(0); //把第一个点删除
                    //上次留下的点重新初始化为无效值
                    pt[0].X = float.NaN;
                    pt[0].Y = float.NaN;
                    pt[1].X = float.NaN;
                    pt[1].Y = float.NaN;
                    CompressLine(ptlst.ToArray(), distance, pt, ref result);
                }

            }

        }

        

        //矢量面的压缩算法：将面分成两条折线再压缩线即可
        //起始点为第一分界点，距起始点最远的点为第二分界点，将面分成两条折线。再对线压缩
        //压缩后再连接，最后再对第二分界点进行判断是留还是删除，根据是第二分界点的前一点和后一点，也是线压缩方法


        /// <summary>
        /// 这是分成两段线的方法。其实还可以直接按一段线来分，尾点增加一次首点即可（唉，光看书了，没想到这个）
        /// </summary>
        /// <param name="polygon"></param>
        /// <param name="distance"></param>
        /// <param name="result"></param>
        private static void CompressPolygon(Coordinate[] polygon, double distance, ref List<Coordinate> result)
        {
            //找出最远的点
            int k = 0;
            double maxDistance = 0;
            for (int i = 1; i < polygon.Length - 1; i++)
            {
                double d = Math.Sqrt((polygon[i].X - polygon[0].X) * (polygon[i].X - polygon[0].X) + (polygon[i].Y - polygon[0].Y) * (polygon[i].Y - polygon[0].Y));
                if (d > maxDistance)
                {
                    maxDistance = d;
                    k = i;
                }
            }

            List<Coordinate> line1 = new List<Coordinate>();
            List<Coordinate> line2 = new List<Coordinate>();

            for (int i = 0; i <= k + 1; i++) // 这里让k加1，也就同时压缩了第二分界点
            {
                line1.Add(polygon[i]);
            }
            for (int i = k; k < polygon.Length; i++)
            {
                line2.Add(polygon[i]);
            }

            List<Coordinate> result1 = new List<Coordinate>();
            List<Coordinate> result2 = new List<Coordinate>();

            //调用压缩线的方法，有两种
            ///第一种
            CompressLine(line1.ToArray(), distance, ref result1);
            CompressLine(line2.ToArray(), distance, ref result2);
            //第二种
            //YaSuoLine(line1.ToArray(), distance, ref result1);
            //YaSuoLine(line2.ToArray(), distance, ref result2);

            result1.RemoveAt(result1.Count - 1); //第一个结果删除最后一个点，因为跟线2是重复的
            //再判断倒数第二个点，即第二分界点在结果1里存不存在，如果存在则线1删除以免重复，线2里面保存就成。如果不存在说明第二分界点不符合条件，线二里面也要删除
            if (result1.Contains(result2[0])) //或者这样判断if(result1[result1.Count -1].X ==result2 [0].X &&result1[result1.Count -1].Y ==result2 [0].Y )                
            {
                result1.RemoveAt(result1.Count - 1);
            }
            else
            {
                result2.RemoveAt(0);
            }

            result.AddRange(result1.ToArray());
            result.AddRange(result2.ToArray());
        }



    }
}
